library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#library(corrplot)
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
source("~/Documents/GitHub/RISKPLOTS/CODE/auxplots.R")
data(cgd)

data <- model.frame(Surv(time,status)~.*.,lung,na.action=NULL)
colnames(data) <-str_replace_all(colnames(data),":","_")
colnames(data) <-str_replace_all(colnames(data),"\\.","_")

data$inst <- NULL
data$`Surv(time, status)` <- NULL
dataLung <- cbind(time=lung$time/365,status=lung$status-1,data)
dataLung <- dataLung[complete.cases(dataLung),]
pander::pander(table(dataLung$status))
0 1
47 121
pander::pander(summary(dataLung$time))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0137 0.4788 0.7356 0.8494 1.14 2.8

Modeling

ml <- BSWiMS.model(Surv(time,status)~1,data=dataLung,NumberofRepeats = 10)

[++++++++++++++++++++++++++]..

sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
  Estimate lower HR upper u.Accuracy
ph_ecog 0.4323 1.194 1.541 1.988 0.6786
sex -0.4591 0.456 0.6319 0.8756 0.6488
pat_karno -0.001767 0.9968 0.9982 0.9997 0.506
ph_karno -3.482e-07 1 1 1 0.5774
Table continues below
  r.Accuracy full.Accuracy u.AUC r.AUC full.AUC
ph_ecog 0.6488 0.6012 0.6012 0.6196 0.5995
sex 0.6786 0.6012 0.6196 0.6012 0.5995
pat_karno 0.7202 0.506 0.5855 0.5 0.5855
ph_karno 0.7202 0.5774 0.57 0.5 0.57
  IDI NRI z.IDI z.NRI Delta.AUC Frequency
ph_ecog 0.04491 0.4048 3.326 2.485 -0.02005 1
sex 0.0285 0.4783 2.758 2.85 -0.00167 1
pat_karno 0.02917 0.3418 2.439 2.243 0.08546 1
ph_karno 0.01425 0.2799 2.221 1.642 0.06998 0.6

Cox Model Performance

Here we evaluate the model using the RRPlot() function.

The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

index <- predict(ml,dataLung)
timeinterval <- 1 # One year

h0 <- sum(dataLung$status & dataLung$time < timeinterval)
h0 <- h0/nrow(subset(dataLung,time > timeinterval | status==1))

rdata <- cbind(dataLung$status,ppoisGzero(index,h0))

rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=dataLung$time,
                     title="Raw Train: Lung Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

As we can see the Observed probability as well as the Time vs. Events are not calibrated.

Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.5463 0.3852 0.2686 0.2686 0.5242
RR 1.214 1.789 60.5 60.5 1.256
RR_LCI 1.014 1.269 0.1278 0.1278 1.028
RR_UCI 1.454 2.523 28638 28638 1.536
SEN 0.314 0.843 1 1 0.6033
SPE 0.8298 0.4894 0.1702 0.1702 0.5957
BACC 0.5719 0.6662 0.5851 0.5851 0.5995
NetBenefit 0.1689 0.5177 0.635 0.635 0.31
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.309 1.086 1.564 0.004168
pander::pander(rrAnalysisTrain$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 121 100.4 144.6 92.43 1.309 1.086 1.564 0.004168
low 83 66.11 102.9 61.73 1.345 1.071 1.667 0.008955
90% 38 26.89 52.16 31.03 1.225 0.8666 1.681 0.208
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
1.03 1.029 0.9971 1.066
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
1.438 1.439 1.424 1.453
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.6516 0.6522 0.5874 0.7095
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.6912 0.5983 0.7841
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.314 0.2327 0.4047
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.8298 0.6919 0.9235
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.5453
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 7.945448 on 1 degrees of freedom, p = 0.004821
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 122 83 95.57 1.654 7.945
class=1 46 38 25.43 6.215 7.945

Cox Calibration

op <- par(no.readonly = TRUE)


calprob <- CoxRiskCalibration(ml,dataLung,"status","time")

( 2.32615 , 781.9918 , 1.468558 , 121 , 140.6966 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
1.555 1.828 2.326

The RRplot() of the calibrated model

h0 <- calprob$h0
timeinterval <- calprob$timeInterval;

rdata <- cbind(dataLung$status,calprob$prob)


rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=dataLung$time,
                     title="Train: Lung",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

Calibrated Train Performance

pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.8523 0.6919 0.531 0.531 0.531
RR 1.214 1.789 60.5 60.5 60.5
RR_LCI 1.014 1.269 0.1278 0.1278 0.1278
RR_UCI 1.454 2.523 28638 28638 28638
SEN 0.314 0.843 1 1 1
SPE 0.8298 0.4894 0.1702 0.1702 0.1702
BACC 0.5719 0.6662 0.5851 0.5851 0.5851
NetBenefit -0.04844 0.2865 0.4575 0.4575 0.4575
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.258 1.044 1.503 0.01424
pander::pander(rrAnalysisTrain$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 121 100.4 144.6 96.17 1.258 1.044 1.503 0.01424
low 83 66.11 102.9 64.23 1.292 1.029 1.602 0.02437
90% 38 26.89 52.16 32.29 1.177 0.8329 1.616 0.2912
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
0.9895 0.9888 0.9577 1.021
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
1.02 1.02 1.015 1.025
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.6519 0.6514 0.5895 0.7068
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.6908 0.5978 0.7838
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.314 0.2327 0.4047
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.8298 0.6919 0.9235
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.8514
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 7.945448 on 1 degrees of freedom, p = 0.004821
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 122 83 95.57 1.654 7.945
class=1 46 38 25.43 6.215 7.945

Crossvalidation

rcv <- randomCV(theData=dataLung,
                theOutcome = Surv(time,status)~1,
                fittingFunction=BSWiMS.model, 
                trainFraction = 0.9,
                repetitions=200,
                classSamplingType = "LOO"
         )
## .[+++].[++].[++].[++].[+++].[++-].[++].[+++].[+++].[++++]10  Tested: 20 Avg. Selected: 3.6 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4271312 
## 
.[++++].[++++].[+++].[++].[++++].[+++].[+++].[+].[++].[++]20  Tested: 40 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.422639 
## 
.[+++].[+++].[++].[++].[+++].[+++].[+++].[+++].[+++].[++]30  Tested: 60 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4261273 
## 
.[++-].[++].[++].[++].[++].[+++].[+++].[++].[++].[+++]40  Tested: 80 Avg. Selected: 3.6 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4419494 
## 
.[+++].[+++].[+++].[+++].[++-].[+++].[+++].[+++].[+++].[+++]50  Tested: 97 Avg. Selected: 3.66 Min Tests: 1 Max Tests: 2 Mean Tests: 1.030928 . MAD: 0.4486325 
## 
.[+++].[+++].[++-].[++].[+++].[+++].[++++].[+++].[+++].[+++]60  Tested: 107 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 2 Mean Tests: 1.121495 . MAD: 0.4708602 
## 
.[++].[+++].[+++].[++].[++].[+++].[++].[+++].[+++].[++-]70  Tested: 117 Avg. Selected: 3.671429 Min Tests: 1 Max Tests: 2 Mean Tests: 1.196581 . MAD: 0.4750843 
## 
.[+].[+++].[+++].[+++].[+++].[++].[+++].[++].[+++].[++]80  Tested: 127 Avg. Selected: 3.65 Min Tests: 1 Max Tests: 2 Mean Tests: 1.259843 . MAD: 0.4704381 
## 
.[++].[+++].[++].[+++].[+++].[+++].[+++].[+++].[+++].[+++]90  Tested: 137 Avg. Selected: 3.666667 Min Tests: 1 Max Tests: 2 Mean Tests: 1.313869 . MAD: 0.4791172 
## 
.[+++].[+++].[+++].[+++].[+++].[+++].[++].[+++].[+++].[++]100  Tested: 147 Avg. Selected: 3.68 Min Tests: 1 Max Tests: 3 Mean Tests: 1.360544 . MAD: 0.4774504 
## 
.[++].[+++].[+++].[++++].[++++].[+++].[+++].[++].[++++].[+++]110  Tested: 157 Avg. Selected: 3.718182 Min Tests: 1 Max Tests: 3 Mean Tests: 1.401274 . MAD: 0.4859387 
## 
.[++].[+++].[++].[+++].[+++].[+++].[+++].[++].[++].[++]120  Tested: 167 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 3 Mean Tests: 1.437126 . MAD: 0.4880085 
## 
.[+++].[++].[+++].[++].[++].[++-].[++++].[++].[++].[+++]130  Tested: 168 Avg. Selected: 3.676923 Min Tests: 1 Max Tests: 3 Mean Tests: 1.547619 . MAD: 0.4865085 
## 
.[+++].[+++].[+++].[+++].[+++].[+++].[+].[+++].[++-].[+++]140  Tested: 168 Avg. Selected: 3.678571 Min Tests: 1 Max Tests: 3 Mean Tests: 1.666667 . MAD: 0.4884544 
## 
.[+++].[+++].[+++].[++-].[+++].[+++].[++].[++].[+++].[++]150  Tested: 168 Avg. Selected: 3.673333 Min Tests: 1 Max Tests: 4 Mean Tests: 1.785714 . MAD: 0.4898956 
## 
.[++++].[+++].[+++].[+++].[++-].[+++].[+++].[++].[+++].[++-]160  Tested: 168 Avg. Selected: 3.68125 Min Tests: 1 Max Tests: 4 Mean Tests: 1.904762 . MAD: 0.4911541 
## 
.[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[++-].[++]170  Tested: 168 Avg. Selected: 3.688235 Min Tests: 1 Max Tests: 4 Mean Tests: 2.02381 . MAD: 0.4897295 
## 
.[+++].[+++].[+++].[++-].[+++].[+++].[+++].[++].[++].[+++]180  Tested: 168 Avg. Selected: 3.688889 Min Tests: 1 Max Tests: 4 Mean Tests: 2.142857 . MAD: 0.4894258 
## 
.[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++]190  Tested: 168 Avg. Selected: 3.705263 Min Tests: 1 Max Tests: 5 Mean Tests: 2.261905 . MAD: 0.4880488 
## 
.[++].[+++].[+++].[++-].[++].[+++].[+++].[+++].[++++].[+++]200  Tested: 168 Avg. Selected: 3.71 Min Tests: 1 Max Tests: 5 Mean Tests: 2.380952 . MAD: 0.4906368 
## 
bbx <- boxplot(rcv$survTestPredictions[,1]~rownames(rcv$survTestPredictions),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(rcv$survTestPredictions[,2]~rownames(rcv$survTestPredictions),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(rcv$survTestPredictions[,4]~rownames(rcv$survTestPredictions),plot=FALSE)$stats[3,],h0)

rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names

rrAnalysisTest <- RRPlot(rdatacv,atRate=c(0.90),
                     timetoEvent=times,
                     title="Test: Lung Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

Calibrated Train Performance

pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.8605 0.6566 0.5676 0.513 0.513
RR 1.186 2.958 3.293 14.58 14.58
RR_LCI 0.972 1.387 1.392 0.03466 0.03466
RR_UCI 1.447 6.309 7.792 6133 6133
SEN 0.1983 0.9587 0.9669 1 1
SPE 0.8936 0.2979 0.2766 0.04255 0.04255
BACC 0.546 0.6283 0.6218 0.5213 0.5213
NetBenefit -0.04053 0.3149 0.4309 0.4382 0.4382
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.241 1.03 1.483 0.01975
pander::pander(rrAnalysisTest$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 121 100.4 144.6 97.5 1.241 1.03 1.483 0.01975
low 98 79.56 119.4 73.58 1.332 1.081 1.623 0.006059
90% 23 14.58 34.51 24.31 0.9463 0.5999 1.42 0.919
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
0.9815 0.9819 0.9504 1.014
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
0.9737 0.974 0.9627 0.9842
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.6027 0.6042 0.5362 0.6737
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.6052 0.5027 0.7077
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.1901 0.1245 0.2714
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.8936 0.769 0.9645
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.8619
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 2.346444 on 1 degrees of freedom, p = 0.125569
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 140 98 103.8 0.3276 2.346
class=1 28 23 17.17 1.981 2.346

Calibrating the test results

rdatacv <- cbind(status,prob,times)
calprob <- CalibrationProbPoissonRisk(rdatacv)
## ( 2.356124 , 781.9918 , 1.487481 , 121 , 143.3503 )
pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.8547 1.005 2.356
timeinterval <- calprob$timeInterval;

rdata <- cbind(status,calprob$prob)


rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90),
                     timetoEvent=times,
                     title="Calibrated Test: Lung",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

### Calibrated Test Performance

pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
Threshold values
  @:0.9 @MAX_BACC @MAX_RR @SPE100 p(0.5)
Thr 0.8618 0.6585 0.5694 0.5147 0.5147
RR 1.186 2.958 3.293 14.58 14.58
RR_LCI 0.972 1.387 1.392 0.03466 0.03466
RR_UCI 1.447 6.309 7.792 6133 6133
SEN 0.1983 0.9587 0.9669 1 1
SPE 0.8936 0.2979 0.2766 0.04255 0.04255
BACC 0.546 0.6283 0.6218 0.5213 0.5213
NetBenefit -0.04266 0.3119 0.4289 0.4362 0.4362
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
O/E Low Upper p.value
1.251 1.038 1.494 0.01672
pander::pander(rrAnalysisTest$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 121 100.4 144.6 96.74 1.251 1.038 1.494 0.01672
low 98 79.56 119.4 73.01 1.342 1.09 1.636 0.00489
90% 23 14.58 34.51 24.12 0.9537 0.6046 1.431 0.9189
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
O/E Mean
mean 50% 2.5% 97.5%
0.9892 0.9892 0.9567 1.022
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
mean 50% 2.5% 97.5%
0.9732 0.9733 0.9627 0.9828
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.6027 0.6017 0.5349 0.6684
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.6052 0.5027 0.7077
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.1901 0.1245 0.2714
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.8936 0.769 0.9645
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90%
0.8632
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 2.346444 on 1 degrees of freedom, p = 0.125569
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 140 98 103.8 0.3276 2.346
class=1 28 23 17.17 1.981 2.346